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Abstract 



We investigate to what extent the specific heat of amorphous sihca can be 
calculated within the harmonic approximation. For this we use molecular 
dynamics computer simulations to calculate, for a simple silica model (the 
BKS potential), the velocity autocorrelation function and hence an effective 
density of states g{i'). We find that the harmonic approximation is valid for 
temperatures below 300K but starts to break down at higher temperatures. 
We show that in order to get a reliable description of the low frequency part 
of g{i')i i.e. where the boson peak is observed, it is essential to use large 
systems for the simulations and small cooling rates with which the samples 
are quenched. We find that the calculated specific heat is, at low temperatures 
(below 50K), about a factor of two smaller than the experimental one. In the 
temperature range 200K <T<Tg, where Tg = 1450K is the glass transition 
temperature, we find a very good agreement between the theoretical specific 
heat and the experimental one. 

PACS numbers: 61.20.Lc, 61.20.Ja, 02.70.Ns, 64.70.Pf 
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I. INTRODUCTION 



It is well known that many low-temperature properties of crystals, such as the specific heat 
or the vibrational dynamics of the atoms, can be calculated if the wave-vector dependence 
of the density of states (DOS) is known. At higher temperatures the system will in general 
become anharmonic and thus the validity of the harmonic approximation breaks down. 
For amorphous systems at low temperatures the situation is similar to the one of crystals. 
However, if the temperature is increased not only anharmonic effects have to be taken into 
account but also the relaxation dynamics of the system since the latter will in general even 
take place at temperatures significantly below the glass transition temperature Tg and thus 
make the harmonic approximation invalid. That the low temperature specific heat Cy of 
silica can indeed be calculated reliably from the DOS was demonstrated by Buchenau et 
al. who used the DOS, as determined from neutron scattering, to calculate Cy between 5 
and 20 Kelvin |Q. Very recently Taraskin and Elliott presented the results of a computer 
simulation in which they had calculated the specific heat of silica in a similar temperature 
range 0. They found that the calculated specific heat is smaller than the experimental 
one and conjectured that the discrepancy might be due to the relatively small size of their 
system. 

What so far has not been investigated is to what extent the harmonic approximation 
can be used to calculate the specific heat also at higher temperatures, i.e. in the range 
100 K < T < Tg, where Tg = 1450 K is the glass transition temperature. In the present 
work we therefore present the results of a calculation of the specific heat in this temperature 
range. For this we use the DOS obtained from a large scale molecular dynamics computer 
simulation and compare the so obtained results with the experimental values. 

The rest of the paper is organized as follows: In the next section we give the details of 
the model used and of the molecular dynamics simulation. In Sec. ^TT|we present the results 
and summarize these in the last section. 



In this section we give the details of the simulations and describe how the glass samples 
have been generated. 

The silica model we are using in this simulation is the one proposed by van Beest et 
al. 0. In this model the interactions (pivij) between two ions i,j that are separated by a 
distance rij is given by 



The values of the partial charges g^, measured in units of e, as well as the values of the 
constants Aij, Bij, and Cij can be found in Refs. One of the remarkable features of this 
potential is that it contains only two-hodj interactions, thus making it very appealing from a 
computational point of view, since the evaluation of computationally demanding three-body 
forces is avoided. Despite its relative simplicity, previous investigations have shown that this 
model gives a realistic description of the static properties of silica glass |^,|] (temperature 
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dependence of the density, structure factor, bond-angle distribution functions), and also of 
the dynamical properties of silica melts |^,|^ (diffusion constants, viscosity, intermediate 
scattering functions). 

For the present simulation we use 8016 ions, in a rigid cubic box of size 48.37A. This 
corresponds to a density of 2.37 g/cm'^ which is close to the experimental value for amorphous 
silica and a pressure which is a bit higher than the normal pressure. The reason for choosing 
a system size that, in the field of simulations of supercooled liquids and glasses, is relatively 
large, is that the dynamics of strong glass formers [§] and also the DOS (see below) shows 
strong finite size effects which have to be avoided. 

The equations of motion have been integrated in the microcanonical ensemble with the 
velocity form of the Verlet algorithm using a time step of 1.6 fs. To generate the glass 
we proceeded as follows: First we equilibrated the system for 4 million time steps at a 
relatively high temperature (2900 K). At this temperature the melt is already quite viscous 
(120 Poise) [^|Tn[ and the typical relaxation times are of the order of 3 ns 0. Subsequently 
the system was coupled to an external heat bath whose temperature was decreased linearly 
in time for 1.0 million time steps to zero Kelvin. This corresponds to a cooling rate of about 
1.8-10^^ K/s. With this cooling rate the system falls out of equilibrium at around 2850 K 
and thus this temperature corresponds to about the value of the fictive temperature of the 
glass. During the cooling procedure we stored periodically copies of the sample. These copies 
were then used as starting configurations for a run in the canonical ensemble with 500,000 
time steps (0.8 ns) in order to anneal the configurations. Afterwards (microcanonical) runs 
were started in which the time Fourier transform of the velocity autocorrelation function 
was measured which was subsequently used to calculate the density of states. In order to 
improve the statistics of the results we averaged over four independent runs. 



III. RESULTS 

In order to calculate the specific heat within the harmonic approximation we need the 
density of states (DOS), g^u). To obtain g{y) one can, e.g., proceed as follows: For any given 
configuration of particles one uses a steepest descent method to determine the location of 
the nearest (metastable) minimum of the potential energy. The DOS can then be calculated 
from the eigenvalues of the (mass-weighted) Hessian matrix at this local minimum. This 
approach was, e.g., used in Ref. 0] to determine how the DOS depends on the cooling rate 
with which the sample was cooled from the high temperature phase to temperatures below 
Tg. Below we will come back to these results. 

A different method to obtain the DOS is to determine the time Fourier transform of the 
velocity autocorrelation function and to make use of the relation |[TT[| 

1 r°° 

9(^^ = j;^T.^^J_^dtexp{z2nut){^^{t)Y,{0)) . (2) 

It should be noted that this approach gives an effective DOS which, however, coincides 
at low temperatures, i.e. when the harmonic approximation is valid, with the real one. 

In Fig. |l| we show this effective DOS for three different temperatures, T = 30 K, T = 
300 K, and T = 1050 K. Note that all of these temperatures are well below the glass 
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transition temperature (Tg=1450K) and thus no relaxational processes take place. From the 
figure we recognize that the main feature of the DOS is a double peak at high frequencies and 
a relatively fiat region at intermediate frequencies, as already discussed in Refs. P,^p!2|. It 
has been shown that the two peaks at high frequencies are due to intra-tetrahedral motions 
of the ions whereas the frequencies of the inter-tetrahedral motions are in the broad region of 
gi^u) at lower frequencies [p!3| , |l4| . From the temperature dependence of the DOS we see that 
the main change in (^(z/) occurs at high frequencies in that, with decreasing temperature, the 
height of the two peaks increases and that their location shifts to higher frequencies, i.e. that 
the corresponding vibrations become faster. (We also note that the curve for T = 1050K 
shows a tail which extends beyond 40THz. We have checked that this feature is not an effect 
of the sampling procedure and thus conclude that it originates from the anharmonicity of 
the system.) At intermediate and low frequencies the temperature dependence of the DOS is 
much weaker. All these observations are in qualitative agreement with the results of Vollmayr 
et al. where a similar dependence was found when the glass transition temperature Tg was 
varied [Q. It should be noted, however, that in that work the dependence of giy) on Tg was 
related to the fact that the structure of the system changes when the cooling rate is varied. 
In contrast to this we investigate here only one cooling rate and the temperature dependence 
of the effective DOS stems only from anharmonic effects and not from a structural change. 

We also note that at very small frequencies the DOS shows a gap. As already discussed 
in Ref. @] this is a finite size effect in that excitations which have a spatial extension larger 
than the size of the simulation box are not present in this system. Some of these excitations 
are, e.g., acoustic phonons with a large wave length. In the Debye theory the density of 
states of these phonons is given by 

goiv) = Su'/ul (3) 

where the Debye frequency is related to q and q, the transverse and longitudinal 
speed of sound, by 



Here is the number of particles in the volume V. By measuring the dispersion relation 
of the transversal and longitudinal acoustic waves at small wave- vectors, we have determined 
the values of q and q [0,0]. We have found that these values, q = 3772m/s and q = 5936m/s, 
agree very well with the experimental values Cf = 3767m/s |jl5| and q = 5970m/s [l^, thus 
demonstrating that with respect to these quantities our model is quite realistic[]. For z/^ we 
obtain 10.65THz which compares well with the experimental value which is 10.40THz ||17|| . 
Therefore we have closed the mentioned gap in our DOS by adding to our measured 5^(1^) 
the Debye-law given by Eq. in the frequency range < z/ < 0.73THz. We will show 



Note that these figures are valid for 300 K. Also it has to be taken into account, that the density of 
our system was fixed at 2.37g/cm'^, which is slightly higher than the experimental value 2.2g/cm^. 
Hence our sound velocities have to be multiplied by ^^2.37/2.2 before they are compared with the 
experimental values. If this is done one obtains q = 3915m/s and q = 6161m/s. 
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below that this modification leads to a significant improvement of the specific heat at low 
temperatures. 

At low temperatures the temperature dependence of the specific heat will be dominated 
by the low frequency part of the spectrum. This frequency range is also of particular interest 
because it includes the so-called boson peak, a dynamical feature whose nature is currently 
a matter of intense debate [p],^ p!6l , |l8| . As mentioned above, at low frequencies the DOS is 
expected to show the mentioned Debye behavior, i.e. glu) oc z/^. In Fig. ^ we therefore 
plot g{v)/v'^ for the three different temperatures (bold lines). We see that the DOS shows 
a large peak at around 1.4THz, i.e. an excess DOS over the Debye value (horizontal dashed 
line), the mentioned boson peak p|,|2|, p!6| , p!8| . Comparing the curves for the three different 
temperatures, we see that in this frequency range the harmonic approximation holds up to 
temperatures of 300K, since only the curve for T = 1050K (dotted line) differs significantly 
from the ones at the two lower temperatures. Note that the curves shown are the one without 
the mentioned Debye correction. In order to see the effect of this correction we include the 
DOS for T = 30K when the correction is taken into account (dash-dotted curve). We see 
that, as expected by construction, this corrected curve goes for small frequencies to the 
Debye value and joins the measured DOS at 0.73THz. 

Also included in the figure are two curves from the work of VoUmayr et al. . These 
were obtained by cooling a sample of = 1002 ions with two different cooling rates 7, 
namely 7 = 7.1 ■ lO^^K/s and 7 = 4.4 ■ lO^^K/s, to zero temperature and calculating the 
DOS from the dynamical matrix. From these two curves we recognize that the height of the 
boson-peak depends on the cooling rate, as already discussed in Ref. |Q, and that the low 
frequency side of the peak of the two curves is at a higher frequency than the one of the 
curves for 8016 ions. One might be tempted to add the mentioned Debye correction also to 
the curves for the smaller system, but a look at the magnitude of these corrections shows 
that they are not sufficient to bring the curves for A^ = 1002 in coincidence with the ones 
of A^ = 8016. Thus we conclude that the excitations giving rise to the boson peak do have 
a significant spatial extension and thus it is necessary to use large system sizes if finite size 
effects have to be avoided. 

From the DOS it is easy to calculate the temperature dependence of Cy, the specific 
heat at constant volume. It is given by 

fc^J ^ "'0 (expi^hu/kBT) - 1) 

At low temperatures the specific heat is expected to show a Debye law, i.e. Cy oc T^, 
and therefore it is useful to plot Cy/T^ vs. T, which is done in Fig. ^. The bold dashed 
curve is the specific heat as obtained from Eq. (^ using the measured DOS at T = 30K and 
the solid bold curve is Cy from this DOS and the Debye corrections. The horizontal dashed 
line is the Debye law 

CyiT) = ^ , (6) 



e 
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with the Debye temperature = huu/kB, which for our model is O^ = 511K. We 
see that the difference between the uncorrected and corrected curve becomes noticeable for 
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temperatures below 20K in that the former decreases towards zero for T ^ whereas the 
latter approaches the Debye value given by Eq. (P). 

Also included in the figure is experimental data by Sosman and by Zeller and Pohl |]19[ . 



We see that although the corrected curve looks qualitatively similar to the experimental 
data, the location of the peak is overestimated by about 5 K. Furthermore also the height 
of the peak, as well as the other parts of the theoretical curve, underestimates the measured 
values by about a factor of two, a result which is in agreement with findings of Taraskin 
and Elliott p|. This discrepancy is not due to the presence of two-level systems, since these 
are relevant only at lower temperatures . Also shown are two curves that stem from the 
DOS calculated by Vollmayr et al. (see Fig. |^). A comparison of these two curves with 
the uncorrected one from the present simulation shows that there are at least two possible 
reasons for the observed discrepancy between the theoretical and experimental curves for 
the specific heat. First of all we see from the curves for N = 1002 ions that there is a 
small but noticable cooling rate dependence of Cy(T) in that the specific heat increases 
with decreasing cooling rate. Thus it can be expected that if our silica model would be 
cooled with a typical experimental cooling rate [0(lK/s)], the theoretical curve would be 
significantly higher than the one shown in the figure. Secondly we see that the theoretical 
curves also show a substantial system size dependence in that the one for the larger system 
is, at low temperatures, significantly above the one for the smaller system. (The fact that 
the curve for the larger system has a cooling rate which is by a factor of two smaller than 
the one for the smaller cooling rate of the smaller system is not sufficient to explain this 
discrepancy between the curves for the two system sizes.) Hence we conclude that the 
anomalous behavior of the specific heat is also considerably affected by finite size effects. 

From the figure we also see that with increasing temperature the relative difference 
between the theoretical and the experimental curves decreases. Therefore it is interesting 
to compare these curves also at higher temperatures, which is done in Fig. ^. Here we 
show Cy(T), as calculated from our DOS at 30K (solid line), and Cp{T), the specific heat 
at constant pressure, as determined by various experiments |]T9| , pT| (lines with symbols). 
We recognize that for temperatures below the glass transition temperature Tg = 1450K 
the agreement between our calculation and the experimental data is very good in that the 
difference is smaller than 0.05 J/gK (see inset). At Tg the theoretical curve starts to deviate 
from the experimental data since at this temperature the real system starts to flow and thus 
the translational degrees of freedom associated with this motion contribute to the specific 
heat. Since the harmonic model does not include this type of motion the theoretical curve 
does not show any special temperature dependence at Tg and with increasing temperature 
approaches the classical value of Dulong and Petit, 1.236J/gK. 

Since the experimental data shown in the figure stem from measurements at constant 
pressure, one has to check whether the discrepancy between the experimental and theoretical 
curves is due to the difference between Cy and Cp. This difference can be expressed by the 
thermodynamic relation 

Cp-Cv = TV— , (7) 

where a is the thermal expansion coefficient and kt is the isothermal compressibility. Using 
the experimental values a = 5.5 ■ 10~^K~^, valid in the temperature range 293K < T < 
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593K [|22|, and kt = 2.79 ■ 10~^GPa~^ at T = 1173K p3[, one obtains an estimate for 
Cp — Cv ~ 5 ■ 10~^J/gK. This is about one decade less than the observed difference between 
our theoretical values and the experimental data and hence we conclude that the observed 
(small) discrepancy is due to some inadequacy of the model. (We also mention that from the 
temperature dependence of the static structure factor at small wave- vectors we have shown 
that the compressibility of our system is in reasonable good agreement with experimental 
values, therefore giving support to the above estimate P,p^. Finally we mention that the 
observed discrepancy also does not disappear if instead of the DOS for 30 K we use the 
one at 300 K or at 1050 K. Hence the difference between the theory and the experiment is 
indeed due to the used model for silica. 



IV. CONCLUSIONS 

We have presented results from a molecular dynamics computer simulation of amorphous 
silica in which we used the velocity autocorrelation function to calculate an effective density 
of states giy). We find that for temperatures less than around 300 K g^p) is essentially 
independent of T, i.e. that the harmonic approximation is valid, whereas significant devia- 
tions are observed for T = 1050 K. Since this last temperature is still significantly below the 
experimental galss transition temperature (Tg=1450 K), and even further below the glass 
transition temperature of this simulation (T^^sim ~ 2850K) we do not expect that the system 
relaxes on the time scale of the simulation. Hence the temperature dependence of g{v) is 
due to the anharmonic nature of the local potential wells in which the ions are sitting. 

We demonstrate that for frequencies in the vicinity of the boson-peak, i.e. around 1- 
2THz, the DOS shows a quite strong dependence on the cooling rate, and more important on 
the system size. Hence it can be concluded that the excitations giving rise to the boson-peak 
have a relatively large spatial extention so that even our simulation box, which measures 
about 48A, is not able to include all of them. Because of this limited system size the low- 
temperature specific heat we obtain from our simulation is about a factor of two smaller 
than the one found in experiments. For higher temperatures we find, however, that our 
theoretical curves agree very well with the experimental data, thus showing that this type 
of calculation is surprisingly reliable for temperatures below Tg. 
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FIGURES 
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FIG. 2. g{u)/v'^ for T =30K (solid line), T = 300K (dotted line) and T = 1050K (dashed line). 
The dashed-dotted curve is the T = 30K curve including the Debye corrections. The two curves 
with symbols are from the simulation of Vollmayr et al. for a = 1002 ion system. The 
horizontal dashed line is the Debye value. 
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FIG. 3. Low temperature dependence of the specific heat. The bold dashed curve is Cy/T^ as 
calculated from the DOS at 30K and the bold curve is obtained if the Debye corrections are taken 
into account. The curves with symbols are the experimental results by Sosman (filled circles) and 
by Zeller and Pohl (open squares) |jT^. The dashed-dotted lines are the ones obtained from the 
DOS of Vollmayr et al. for two different cooling rates B. The horizontal line is the Debye law. 



1.5 



o> 1.2 



Sim. from g(v) 



• — • Sosman 
G^-o Zeller/Pohl 
Q — Q Richet et al. 



— ^ ^ — □ B ■ 




900 1200 1500 
T[K] 



800 1200 1600 2000 

T [K] 



FIG. 4. Temperature dependence of the specific heat at constant volume as predicted by the 
harmonic approximation (solid line). The symbols are experimental data for the specific heat at 
constant pressure by Sosman, by Zeller and Pohl, and by Richet et al. |19,^|. Inset: Difference 
between the data of Richet et al. and Sosman and our data. 
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